Integration of Small RNA and Degradome Sequencing Reveals the Regulatory Network of Al-Induced Programmed Cell Death in Peanut

Peanut is one of the most important oil crops in the world. In China, the peanut is highly produced in its southern part, in which the arable land is dominated by acid soil. At present, miRNAs have been identified in stress response, but their roles and mechanisms are not clear, and no miRNA studies have been found related to aluminum (Al)-induced programmed cell death (PCD). In the present study, transcriptomics, sRNAs, and degradome analysis in the root tips of two peanut cultivars ZH2 (Al-sensitive, S) and 99-1507 (Al-tolerant, T) were carried out. Here, we generated a comprehensive resource focused on identifying key regulatory miRNA-target circuits that regulate PCD under Al stress. Through deep sequencing, 2284 miRNAs were identified and 147 miRNAs were differentially expressed under Al stress. Furthermore, 19237 target genes of 749 miRNAs were validated by degradome sequencing. GO and KEGG analyses of differential miRNA targets showed that the pathways of synthesis and degradation of ketone bodies, citrate cycle (TCA cycle), and peroxisome were responded to Al stress. The combined analysis of the degradome data sets revealed 89 miRNA-mRNA interactions that may regulate PCD under Al stress. Ubiquitination may be involved in Al-induced PCD in peanut. The regulatory networks were constructed based on the differentially expressed miRNAs and their targets related to PCD. Our results will provide a useful platform to research on PCD induced by Al and new insights into the genetic engineering for plant stress response.


Introduction
Aluminum (Al) toxicity is a worldwide agronomic problem that limits crop production in acidic soils. The response of plant roots to Al is mainly caused by Al exclusion mechanisms and Al tolerance mechanisms [1]. The mechanism of Al exclusion mainly includes the secretion of organic acids (citrate, malate, and oxalate) [2], the increase of rhizosphere pH [3], and the binding of the cell wall to Al [4]. The mechanisms of Al tolerance mainly include the regionalization of Al, the inactivation of toxicity after the chelation of Al with organic acids in cells, in order to reduce the toxicity of Al to plants and improve plant resistance [5]. It was found that Al stress could induce the changes in antioxidant enzyme activities in plants, including superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), ascorbate peroxidase (APX), etc. [6].
Al can induce plant cells to produce reactive oxygen species (ROS) and activate some antioxidant enzymes. A high concentration of ROS resulted in membrane lipid peroxidation sequencing analyses in peanut. The combination of sRNA, degradome, and transcriptome has produced a comprehensive resource, focusing on identifying the key regulatory pathways of miRNAs that could induce PCD under Al stress and provide useful information for the improvement of peanut plants.

Sequencing and Identification of Known and Novel miRNAs
ZH2 (Al-sensitive, S) and 99-1507(Al-tolerant, T) genotypes were used to generate 26 small RNA libraries representing three controls without Al stress at different time points (0, 8, and 24 h-0) and two Al treatments at different time points (8 and 24 h) for high-throughput sequencing ( Table S1). The bioinformatics analysis of our data was performed using the development pipeline ( Figure S1). After removing low-quality sequences, 18 to 25 nt long sequences were obtained. The size distributions of the total and unique sRNAs were summarized in Figure S2. There was no significant difference among the 26 libraries. The displayed length of these peanut sRNAs ranged from 18 to 25 nt, and the sRNA with 24 nt showed the highest redundancy size category among all of the 26 libraries ( Figures S2 and 1). It revealed that 65.67% of all the miRNAs were 24 nt in length ( Figure 1). ZH2 (Al-sensitive, S) and 99-1507(Al-tolerant, T) genotypes were used to generate 26 small RNA libraries representing three controls without Al stress at different time points (0, 8, and 24 h-0) and two Al treatments at different time points (8 and 24 h) for highthroughput sequencing (Table S1). The bioinformatics analysis of our data was performed using the development pipeline ( Figure S1). After removing low-quality sequences, 18 to 25 nt long sequences were obtained. The size distributions of the total and unique sRNAs were summarized in Figure S2. There was no significant difference among the 26 libraries. The displayed length of these peanut sRNAs ranged from 18 to 25 nt, and the sRNA with 24 nt showed the highest redundancy size category among all of the 26 libraries ( Figures  S2 and 1). It revealed that 65.67% of all the miRNAs were 24 nt in length ( Figure 1).
We used the BLAST algorithm to query all of the clean and unique sRNA read sequences and searched for related non-coding RNAs stored in Rfam and NCBI GenBank databases then removed rRNA, tRNA, snRNA, and snoRNA sequences. In order to predict new miRNAs, the remaining unannotated sRNA unique reads that do not match with any transcriptome library were searched, and the secondary structure prediction was performed. All of the sRNAs loci that could be folded into secondary structure were considered as potential new miRNA candidates. A total of 2284 miRNAs were identified, including 1670 novel (gp4), 83 known (gp1a), 185 potential novel miRNAs (predicted but not mapped in the plant genome, gp1b, 2a), and 346 unknown miRNAs (mapped on the plant genome but not in the miRNA database, gp2b, 3; Figure 2 and Table S2). These miRNAs belong to 49 recorded miRNA families ( Figure S3).  We used the BLAST algorithm to query all of the clean and unique sRNA read sequences and searched for related non-coding RNAs stored in Rfam and NCBI GenBank databases then removed rRNA, tRNA, snRNA, and snoRNA sequences. In order to predict new miRNAs, the remaining unannotated sRNA unique reads that do not match with any transcriptome library were searched, and the secondary structure prediction was performed. All of the sRNAs loci that could be folded into secondary structure were considered as potential new miRNA candidates. A total of 2284 miRNAs were identified, including 1670 novel (gp4), 83 known (gp1a), 185 potential novel miRNAs (predicted but not mapped in the plant genome, gp1b, 2a), and 346 unknown miRNAs (mapped on the plant genome but not in the miRNA database, gp2b, 3; Figure 2 and Table S2). These miRNAs belong to 49 recorded miRNA families ( Figure S3). The four indicated groups of miRNAs identified from an incompatible interaction. The miRNA counts, the counts of miRNAs from miRBase. Expression level, low indicates <10, middle indicates >10, but less than average, high indicates over average. Hairpin length, pre-miRNA length.
In addition to the validation of differentially expressed miRNAs, we also examined the distribution of 147 differential miRNAs between the control and treatments ( Figure  4c). Compared with the control, the largest number (92 miRNAs) was upregulated miR-NAs in ZH2 after 24 h of Al treatment, while 35 miRNAs were significantly downregulated in 99-1507 after 24 h of Al treatment (Figure 4c). It was found that the same miRNA could show different expression patterns up or downregulation in different comparisons. For example, PC-3p-46954_232 was upregulated at 8 h in ZH2 but downregulated at 24 h in ZH2, PC-5p-85555_103 was downregulated at 8 h in ZH2 but was upregulated at 24 h in ZH2 (Figure 4b). These results indicated that the expression levels of miRNAs may depend on the duration of Al treatment.
In addition to the validation of differentially expressed miRNAs, we also examined the distribution of 147 differential miRNAs between the control and treatments (Figure 4c

Target Prediction of the Known and Novel miRNAs by Degradome Sequencing
To research on the functions of miRNA, we first determined the targets of miRNA. Through degradome sequencing, 7.97 million raw reads representing 3.39 million unique reads, were generated from the six mixed degradome pools (Table S4). After removing the reads lacking the adaptor and the low-quality sequences, 5.66 million (71.0% of all reads) sequences were successfully mapped to peanut, containing 2.16 million unique reads (Table S4).
Through degradome analysis, a total of 19237 targets were identified for 749 miRNAs (Table S5). Generally, many miRNAs targeted multiple genes and most of them were highly conserved. Different members of miRNA families cleaved the same target genes or different members of the same gene family. For instance, 68 targets were identified to be cleaved by gma-miR172g_R+1_1ss9AT, and 67 targets were identified to be cleaved by mtr-miR1507-5p_R-1_2ss5TA18AG, while eight miRNAs were detected to cleave only one transcript target (XM_016096779.2), which was annotated to a protein related to the heat shock. In addition, degradome analysis showed that 8192 genes were the targets of 147 differentially expressed miRNAs in peanut (Table S6). Through the expression analysis of target genes in response to Al stress in peanut, a large number of members of the WRKY, ABC transporters, NAC, MYB, and TCP families were ranked in the context of Al response, indicating their important roles for peanut in combating Al stress (Table S6). Generally, many miRNAs targeted PCD-related genes and stress-related genes. For instance, metacaspases (MCs) observed in the study were targeted by both gma-miR396a-5p_1ss21GA and mtr-MIR2592bj-p3_2ss11TC18AT, and miR10420 members targeted SODs, PODs, RBOHs, and APXs, etc. (Table S7).

Annotation and Enrichment Analysis of Targets for miRNAs under Al Stress
To research the regulatory functions of miRNAs in response to Al stress in peanut, gene ontology (GO) functional classification and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment were conducted to further understand the functions ( Figure 5). Out of the 749 genes, differential miRNAs targeted 8215 identified target genes, and the GO terms were assigned to 6724 genes. For these 6724 genes, a total of 4387 GO terms were obtained. It was seen that the GO terms for differential miRNAs were uniformly assigned to each of the biological processes (2540), molecular function (553), and cellular component (1294) categories.
Second, based on the KEGG analysis, the most enriched pathways in response to Al stress were those related to the synthesis and degradation of ketone bodies (ko00072), tyrosine metabolism (ko00350), tryptophan metabolism (ko00380), butanoate metabolism (ko00650), and carbon metabolism (ko00710) (Figure 5b). Previous studies have inferred that miR159, miR164, and miR319 were related to Al-induced PCD [26]. In this study, their target genes are mainly enriched in enzyme activator activity, gluconeogenesis, GTPase activator activity, monosaccharide biosynthetic process, nucleoside-triphosphatase regulator activity, carbohydrate biosynthetic process, etc.

Regulatory Analysis of miRNAs Related to PCD and Target Genes
To further investigate the association of Al-induced PCD responsive miRNAs with their targets, we used the Cytoscape platform to analyze the differential miRNAs and their target genes that may be related to PCD, including endoplasmic reticulum stress, oxidative stress, and mitochondrial pathway (Table S7, Figure 6). Eighty-nine pairs of miRNA-mRNA were identified including 69 target genes and 45 miRNAs. The target genes included metacaspases, cytochrome c, cysteine protease, etc. The results indicated that many miRNAs targeted PODs and PDILs, which were involved in oxidative stress and might participate in Al-induced PCD by regulating ROS. Among them, miR10420 and miR6300 jointly regulated PODs. The miR10420 and miR6300 targeted a variety of mRNAs that may be involved in Al-induced PCD, including SODs and PODs, etc. In addition, miR2592 and miR396 jointly regulated the metacaspases (MCs), which are crucial executioners in Al-induced PCD (Table S7, Figure 6).

Regulatory Analysis of miRNAs Related to PCD and Target Genes
To further investigate the association of Al-induced PCD responsive miRNAs their targets, we used the Cytoscape platform to analyze the differential miRNAs and target genes that may be related to PCD, including endoplasmic reticulum stress, o tive stress, and mitochondrial pathway (Table S7, Figure 6). Eighty-nine pairs of miR mRNA were identified including 69 target genes and 45 miRNAs. The target gen cluded metacaspases, cytochrome c, cysteine protease, etc. The results indicated many miRNAs targeted PODs and PDILs, which were involved in oxidative stres might participate in Al-induced PCD by regulating ROS. Among them, miR10420 miR6300 jointly regulated PODs. The miR10420 and miR6300 targeted a varie mRNAs that may be involved in Al-induced PCD, including SODs and PODs, etc. I dition, miR2592 and miR396 jointly regulated the metacaspases (MCs), which are c executioners in Al-induced PCD (Table S7, Figure 6).  Table S7.

Validation of miRNA and Its Target Gene by RT-qPCR
To confirm the reliability of RNA-Seq data, a validation experiment on random lected genes was determined using qRT-PCR (Figure 7). Expression patterns o Figure 6. The miRNA-mediated regulatory networks responses to Al stress in peanut. Red circles represent the differentially expressed miRNAs, and blue circles represent the target genes of differentially expressed miRNAs. Abbreviations for target genes are listed in Table S7.

Validation of miRNA and Its Target Gene by RT-qPCR
To confirm the reliability of RNA-Seq data, a validation experiment on randomly selected genes was determined using qRT-PCR (Figure 7). Expression patterns of six miRNAs and six mRNAs under Al stress were consistent with those obtained by the deepsequencing, which confirmed that our high-throughput sequencing data and subsequent interpretations were reliable. miRNAs and six mRNAs under Al stress were consistent with those obtained by the deepsequencing, which confirmed that our high-throughput sequencing data and subsequent interpretations were reliable.     To standardize RNA quantity for evaluating relative expression levels, U6 and actin were used as the internal reference genes of miRNAs and their targets, respectively. The expression level at 0 h in each sample was used as reference state, which was set to 1, and fold change values were shown here.

Discussion
Recently, some miRNAs from peanut were identified by computational and direct cloning approaches [27]. Ma et al. identified 1082 miRNAs in developing peanut seeds by high-throughput sequencing [28], but the characteristics and functions of most peanut miRNAs are still unknown. Here, we identified 2284 miRNAs in peanut root tips (Table  S2). The sRNA length distribution patterns peaked at 24 nt in our study (Figures 1 and  S2), which was consistent with previous results for most angiosperms [29]. We found that the sRNA length distribution pattern in peanut was similar to the other plants.
In this study, we identified 147 miRNAs differentially regulated by Al exposure (Table S3). Most of the miRNA members of the same family showed similar expression profiles, for instance, six miR160 and eight miR408 family members were significantly upregulated by Al exposure (Figure 4a,b, Table S3). It has been reported that miR160 could also be upregulated by Al stress in rice [30]. On the contrary, miR160 has been found to be downregulated in barrel medic in response to Al stress [14]. In addition, the miR408 was downregulated by Al exposure in rice [11]. The differential expression patterns of miR160 and miR408 may depend on the different Al treatment conditions and plant species. We also discovered that some previously known miRNAs were regulated by Al stress. Fiftynine miR6300, 10 miR2592, and six miR530 family members were found to be significantly upregulated by Al stress (Figure 4a,b, Table S3). In sugarcane, miR2592 was likely to be the regulator of internode elongation genes [31]. In plants, miR530 was found to be responsive to abiotic stresses, such as salinity [31]. Our findings implied that the known miRNAs were most likely involved in cross adaptation of plants to abiotic stresses. In our results, miR6300 targeted antioxidant enzymes, including SODs, PODs, and CATs. In addition, miR2592 targeted MCs, which might activate caspases and induce PCD. Moreover, miR530 targeted CSE, which might affect oxidative stress to participate in Al-induced PCD. Their exact functions need to be verified in the future studies. In the meanwhile, different miRNA members within the same family were present at different expression levels after Al exposure. For example, mtr-miR390-R+1 was downregulated by Al exposure, while mtr-miR390_1ss21CA was significantly upregulated (Figure 4a,b, Table S3). Similar events also occurred with miRNA family members in other plants. In radish, miR172a was significantly upregulated, while miR172c was downregulated by Cd exposure [32]. Therefore, a complex mechanism may regulate the expression profiles of members of the same family.
Recent studies have shown that many miRNAs regulate the expression of target genes through transcript cleavage or translation inhibition [33]. A transcriptome-wide analysis of the degradome was performed, and numerous target transcripts for the known and novel miRNAs were determined. This study presented the first transcriptome-based analysis of miRNA targets responsive to Al stress in peanut by degradome sequencing To standardize RNA quantity for evaluating relative expression levels, U6 and actin were used as the internal reference genes of miRNAs and their targets, respectively. The expression level at 0 h in each sample was used as reference state, which was set to 1, and fold change values were shown here.

Discussion
Recently, some miRNAs from peanut were identified by computational and direct cloning approaches [27]. Ma et al. identified 1082 miRNAs in developing peanut seeds by high-throughput sequencing [28], but the characteristics and functions of most peanut miRNAs are still unknown. Here, we identified 2284 miRNAs in peanut root tips (Table S2). The sRNA length distribution patterns peaked at 24 nt in our study (Figures 1 and S2), which was consistent with previous results for most angiosperms [29]. We found that the sRNA length distribution pattern in peanut was similar to the other plants.
In this study, we identified 147 miRNAs differentially regulated by Al exposure (Table S3). Most of the miRNA members of the same family showed similar expression profiles, for instance, six miR160 and eight miR408 family members were significantly upregulated by Al exposure (Figure 4a,b, Table S3). It has been reported that miR160 could also be upregulated by Al stress in rice [30]. On the contrary, miR160 has been found to be downregulated in barrel medic in response to Al stress [14]. In addition, the miR408 was downregulated by Al exposure in rice [11]. The differential expression patterns of miR160 and miR408 may depend on the different Al treatment conditions and plant species. We also discovered that some previously known miRNAs were regulated by Al stress. Fifty-nine miR6300, 10 miR2592, and six miR530 family members were found to be significantly upregulated by Al stress (Figure 4a,b, Table S3). In sugarcane, miR2592 was likely to be the regulator of internode elongation genes [31]. In plants, miR530 was found to be responsive to abiotic stresses, such as salinity [31]. Our findings implied that the known miRNAs were most likely involved in cross adaptation of plants to abiotic stresses. In our results, miR6300 targeted antioxidant enzymes, including SODs, PODs, and CATs. In addition, miR2592 targeted MCs, which might activate caspases and induce PCD. Moreover, miR530 targeted CSE, which might affect oxidative stress to participate in Al-induced PCD. Their exact functions need to be verified in the future studies. In the meanwhile, different miRNA members within the same family were present at different expression levels after Al exposure. For example, mtr-miR390-R+1 was downregulated by Al exposure, while mtr-miR390_1ss21CA was significantly upregulated (Figure 4a,b, Table S3). Similar events also occurred with miRNA family members in other plants. In radish, miR172a was significantly upregulated, while miR172c was downregulated by Cd exposure [32]. Therefore, a complex mechanism may regulate the expression profiles of members of the same family.
Recent studies have shown that many miRNAs regulate the expression of target genes through transcript cleavage or translation inhibition [33]. A transcriptome-wide analysis of the degradome was performed, and numerous target transcripts for the known and novel miRNAs were determined. This study presented the first transcriptome-based analysis of miRNA targets responsive to Al stress in peanut by degradome sequencing analysis. For the miRNAs, 19237 targets were identified by degradome sequencing (Table S5). An integrated analysis of miRNA expression profiles and their targets could help in identifying the functional miRNA-target modules involved in regulating specific biological processes, such as Al stress and PCD [34]. Consistent with some previous studies, several identified miRNA targets in peanut, belong to a variety of Al-responsive genes, including WRKY factors, ABC transporter, NAC-domain, and F-box proteins (Table S6) [35,36]. It has been previously shown that miR394 targets mRNA of the F-box gene [37]. In this study, ahy-miR394 targeted WRKY transcription factors, ABC transporter, and F-box proteins. WRKY47 has been demonstrated to be required for Al tolerance and root growth via the regulation of cell wall modification genes [38]. AtSTAR1 was required for Al tolerance in Arabidopsis, probably functioning as a bacterial-type ABC transporter by forming a complex with ALS3 [39]. Fang et al. found the degradation of STOP1 mediated by the F-box proteins RAH1 and RAE1 balanced Al tolerance and plant growth [40]. In addition, ahy-miR3508 targeted ABC transporter and NAC-domain in this study. Moreover, it was reported that the positive effect of Al on the growth of rice is mediated by NACs that respond to hormones [41]. The miRNAs involved in the control of gene expression related to PCD were analyzed in our study. Metacaspases (MCs) are cysteine-dependent proteases, which play essential roles in PCD [42]. The mtr-MIR2592bj-p3_2ss11TC18AT and gma-miR396a-5p_1ss21GAg targeted metacaspases, and several miRNAs including miR10420 and miR2592 targeted RBOHs (Table S7). Overexpression of alternative oxidases (AOXs) alleviated mitochondria-dependent PCD induced by Al phytotoxicity in Arabidopsis [43]. In our results, ptc-MIR171m-p5_1ss17CT and PC-3p-24641_483 also targeted AOXs (Table S7).
Recently, several reports have provided evidence that endoplasmic reticulum stress [44], oxidative stress [45], and mitochondrial pathway [46] are associated with PCD. The presence of many signal pathways shows that plants respond to Al exposure through complex mechanisms. By screening the functional annotations of target genes for differential miRNA, the target genes whose functional annotation related to PCD and their corresponding miR-NAs were identified. A total of 89 miRNA-mRNA pairs were found (Table S7). The results showed that miRNA targeted genes related to endoplasmic reticulum stress, oxidative stress, and mitochondrial pathway. Other target genes were involved in Al-induced PCD through ubiquitination. It has been shown that the ubiquitin-26s proteasome system (UPS) is an important protein degradation system in cells. Under the joint of ubiquitin-activating enzyme (E1), ubiquitin-conjugating enzyme (E2), and ubiquitin ligase (E3), the system binds the tag protein ubiquitin composed of 76 amino acids to the substrate protein. The labeled protein is degraded by proteasome or regulated by signal pathway. Different forms of ubiquitination modification of substrate proteins can regulate a variety of cell activities, such as signal transduction, gene transcription, and PCD [47]. Several studies showed that UPS plays a key and complex role in the regulation of PCD [48]. Among them, the RING finger E3 family (RING finger) has the closest relationship with PCD [49]. In our results, ath-miR159b-3p_1ss7TA and csi-miR167c-3p_1ss10TA might target the ubiquitinconjugating enzyme E2 in response to Al-induced PCD in peanut. The miR159 could respond to Al stress [50] and inhibited GAMYBs transcription factor by inhibiting GAregulated PCD in tapetum cells [51]. The miR167 negatively regulated ARF, which encoded auxin response transcription factors [52]. In B. napus, BnNRAMP1b was a target of miR167 under Cd stress [53]. The miR482 may be involved in Al-induced PCD by targeting E3 ubiquitin-protein ligase SINAT5 (Table S7), it was found to be involved in the defense regulation of plant pathogens [54]. The miR10420 may be involved in Al-induced PCD by targeting 26S protease regulatory subunit and RING finger protein (Table S7). There are few studies on miR10420.
From our findings, miRNA may be involved in Al stress by targeting ubiquitinationrelated genes and then activating caspase in coordination with the mitochondrial pathway to induce PCD. However, some specific mechanisms are still controversial, which will be a hot issue in the future. In addition, most of the previous studies of miRNA in Al stress focused on humans and other animals [55], but few are available on plants. The inhibition of proteasome function by gene silencing could also cause PCD in tobacco cells [56]. The ubiquitin-26s proteasome system of ubR48 (ubiquitin variant, the 48th lysine of the translated peptide chain is replaced by arginine) in transgenic Arabidopsis thaliana was inhibited, which induced PCD [57]. There is an important role of UPS in PCD regulation. However, the regulatory process and its components in the control of PCD in plants should be studied in the future.

Plant Materials and Al Stress Treatment
Two different peanut (Arachis hypogaea L.) varieties 99-1507 (T, Al-tolerant) and ZH2 (S, Al-sensitive) were selected as plant materials. The seeds were first incubated in the moist perlite under dark conditions at 26 ± 1 • C for 3-4 days for germination. After germination, the seed coats were removed, then seedlings with about 2-3 cm roots length were transferred in a modified Hoagland nutrient solution. The growing conditions were set as 26 ± 2 • C, 12 h light/12 h dark and the light intensity was 30-50 µM·m −2 ·s −1 . Seedlings with the third true leaf were pretreated with 0.1 mM CaCl 2 solution (pH 4.2) for 24 h, then treated with 100 µM AlCl 3 solution (containing 0.1 mM CaCl 2 , pH 4.2) for 0, 8, and 24 h. Thereafter, they were treated with 0.1 mM CaCl 2 for the same time (8 and 24 h-0) as the mock. Approximately 1 cm of the root tips were gained for the sample collection. Each treatment was carried out in triple biological replicates except the control 8 and 24 h-0, which were repeated twice. All of the samples were immediately frozen in liquid nitrogen and stored at −80 • C.

The sRNA Sequencing and miRNAs Identification
Twenty six small RNA libraries were constructed using the TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, USA), according to the manufacturer's instructions. Total RNA, which was a qualified sample, was ligated to the RNA 3 and RNA 5 adapters and reversed transcription, followed by PCR to make cDNA constructs of the small RNAs. Finally, the cDNA constructs were purified and the library was validated. Then, single-end sequencing (50 bp) was performed on an Illumina Hiseq2500 at the Hangzhou LC-Bio Co., Ltd. (Hangzhou, China). Thereafter, the dataset was processed using an in-house program, ACGT101-miR (LC Sciences, Houston, TX, USA) to remove adapter dimers, low complexity reads, common RNA families (rRNA, tRNA, snRNA, snoRNA), and RepBase database (http://www.girinst.org/repbase, accessed on 27 December 2021, repetitive sequence database). Subsequently, unique sequences with 18~25 nucleotides in length were mapped to specific species precursors in miRBase 22.0 (http://www.mirbase.org/, accessed on 27 December 2021) by BLAST search to identify known miRNAs and novel 3p-and 5pderived miRNAs.

Analysis of Differentially Expressed miRNAs
To discover the expression profiles of the target genes, 26 independent libraries were constructed from the RNA samples of two varieties from five different Al treatments. For each library, all of the sequences were processed to filter out adaptor and lower quality sequences. Then, all of the valid data were performed for further miRNA identification and prediction analysis. Differential expression of miRNAs based on normalized deepsequencing counts was analyzed by selectively using the Fisher exact test, Chi-squared 2 × 2 test, Chi-squared n × n test, Student t-test or ANOVA based on the experiments design. Then, a rigorous algorithm method was performed to identify the differentially expressed genes between two samples. The significantly differentially expressed genes among all of the different samples were judged by the threshold as follows: p-value < 0.05 and the absolute value of log2 ratio ≥ 1. Differentially expressed miRNAs were identified using UpSetR (v1.4.0) (LC-Bio Co., Ltd., Hangzhou, China) in R. Through mock, the influence of unrelated factors was eliminated. The heatmap of the differentially expressed miRNAs in ZH2 and 99-1507 was constructed using the ggplot2 package (v1.0.12) in R.

Degradome Sequencing, Target Identification, and Analysis
An equal number of three biological replicates of RNA samples at the same time point were mixed together to generate six degradation libraries. Then, 26 sRNA samples and six degradation mixed samples were sent to Hangzhou LC-Bio Co., Ltd. (Hangzhou, China). The extracted sequencing reads generated from degradome sequencing were used to identify potential cleavage targets through the CleaveLand (v4.0). Then, the degradome reads were mapped to the peanut transcriptome data, which was previously tested in our laboratory [24].

Gene Annotation, Functional Enrichment, and Pathway Enrichment Analysis
Expression analysis of all the target genes was performed using transcriptomic data of peanut in response to Al stress [24]. GO functional classification of all the identified targets of differentially expressed miRNAs was performed to uncover the miRNA-gene regulatory network. The KEGG pathway for the target genes was analyzed. Target genes were annotated based on the Gene Ontology database (https://ftp.ncbi.nlm.nih.gov/gene/ DATA/gene2go.gz, accessed on 27 December 2021). Pathway analyses of target genes were performed using The Kyoto Encyclopedia of Gene and Genome (KEGG) database (http://www.genome.jp/kegg). Differential miRNA-mRNA related to Al-induced PCD was presented by functional annotation and the miRNA-mRNA network construction was presented using Cytoscape (v3.4.0). For network construction, stress-responsive miRNAs and genes, whose functional annotation were involved in PCD response, were included.

Verification by qRT-PCR
For validation of miRNAs and their targets expression, high-throughput sequencing technology and qRT-PCR of randomly selected miRNAs and genes were performed. The qRT-PCR analysis was carried out for six miRNAs and six targets at three time points (0, 8, and 24 h). The qRT-PCR analysis was performed with the miRNA 1st Strand cDNA Synthesis Kit (Vazyme Biotech, Nanjing, China) and SYBR Premix Ex TaqTM (TaKaRa, Dalian, China), according to the manufacturer's instructions. RT-qPCR was conducted on a Bio-Rad iQ5 Real-time PCR System (Bio-Rad Laboratories, Hercules, CA, USA). U6 was used as the internal control for miRNAs and actin was used as the internal control for target genes to standardize the RNA quantity for evaluating relative expression levels. All of the primers are listed in Table S8. At least three technical replicates were carried out for each sample and three independent biological replicates were performed.

Conclusions
In summary, we conducted an integrated analysis of the sRNAs and degradome to generate a comprehensive resource focused on identifying the key regulatory miRNA-target for Al-induced PCD. A total number of 2284 miRNAs were identified in peanut. In addition, 147 miRNAs were significantly differentially expressed and identified as Al-responsive miRNAs. We screened the target genes related to endoplasmic reticulum stress, PCD, oxidative stress, etc. Eighty nine pairs of miRNA-mRNA composed of 45 differential miR-NAs and 69 target genes were identified ( Figure 6, Table S7), which may be the key genes related to Al-induced PCD. Based on our small RNA, degradome sequencing, previous transcriptomics analysis, and results of other peanut crops, a putative model of regulatory network associated with Al-induced PCD in peanut was proposed (Figure 9). Al stress triggered H 2 O 2 bursts, which increased calcium concentrations to the point of overload. This triggered ROS production in mitochondria, opened MPTP, released cytochrome c, activated caspases activity, and finally induced the occurrence of PCD. Many miRNAs are involved in this process via their regulation of mRNAs related to PCD. For example, miR10420 targeted their target genes to participate in oxidative stress, endoplasmic reticulum stress, and Cyt c releasing, etc. Through endoplasmic reticulum stress, and the miR2592 targeted MCs are involved in Al-induced PCD. Herein, our results will help in elucidating the miRNA-mediated molecular mechanisms of Al-induced PCD in peanut and in providing insights for plants adaption to Al stress. are involved in this process via their regulation of mRNAs related to PCD. For example, miR10420 targeted their target genes to participate in oxidative stress, endoplasmic reticulum stress, and Cyt c releasing, etc. Through endoplasmic reticulum stress, and the miR2592 targeted MCs are involved in Al-induced PCD. Herein, our results will help in elucidating the miRNA-mediated molecular mechanisms of Al-induced PCD in peanut and in providing insights for plants adaption to Al stress.